WGCNA Results into Functional Enrichment and Biomineralization Toolkit Analysis
Kristen Brown
11/01/2023
# load libraries - notes show the install command needed to install (pre installed)
library(goseq)### Loading required package: BiasedUrn## Loading required package: geneLenDataBase## library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forcats)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.1library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(tidyr)
library(grDevices)
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary(Rmisc)## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.1## Loading required package: plyr## Warning: package 'plyr' was built under R version 4.3.1## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizelibrary(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutatelibrary(tibble)
#library(hrbrthemes)#
library(gridExtra)
library(tidyr)
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric#BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)## Loading required package: grid## ========================================
## ComplexHeatmap version 2.16.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================library(circlize)## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================#BiocManager::install("GSEABase")
library(GSEABase)## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: annotate## Loading required package: AnnotationDbi## Warning: package 'AnnotationDbi' was built under R version 4.3.1## Loading required package: stats4## Loading required package: IRanges## Warning: package 'IRanges' was built under R version 4.3.1## Loading required package: S4Vectors## Warning: package 'S4Vectors' was built under R version 4.3.1##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:plyr':
##
## rename## The following object is masked from 'package:tidyr':
##
## expand## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:geneLenDataBase':
##
## unfactor## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname##
## Attaching package: 'IRanges'## The following object is masked from 'package:plyr':
##
## desc## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: XML## Loading required package: graph##
## Attaching package: 'graph'## The following object is masked from 'package:XML':
##
## addNode## The following object is masked from 'package:circlize':
##
## degree## The following object is masked from 'package:plyr':
##
## joinlibrary(data.table)##
## Attaching package: 'data.table'## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second## The following objects are masked from 'package:reshape2':
##
## dcast, melt## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary(stringr)##
## Attaching package: 'stringr'## The following object is masked from 'package:graph':
##
## boundarylibrary(GenomicRanges)## Warning: package 'GenomicRanges' was built under R version 4.3.1## Loading required package: GenomeInfoDb## Warning: package 'GenomeInfoDb' was built under R version 4.3.1library(rrvgo)## Warning: package 'rrvgo' was built under R version 4.3.1library(rtracklayer)## Warning: package 'rtracklayer' was built under R version 4.3.1geneInfo <- read.csv("../../../output/Slope_Base/geneInfo.csv")
WGCNA_ModuleMembership <- read.csv("../../../output/WGCNA/WGCNA_ModuleMembership.csv") %>% rename(c("X"="gene_id")) #this file was generated from the WGCNA analyses and contains the modules of interest
geneInfo <- geneInfo %>% left_join(WGCNA_ModuleMembership %>% dplyr::select(gene_id, moduleColor), by = "gene_id")Get gene length information.
#import file
gff <- rtracklayer::import("../../../data/Pocillopora_acuta_HIv2.genes_fixed.gff3") #if this doesn't work, restart R and try again
transcripts <- subset(gff, type == "transcript") #keep only transcripts , not CDS or exons
transcript_lengths <- width(transcripts) #isolate length of each gene
seqnames <- transcripts$ID #extract list of gene id
lengths <- cbind(seqnames, transcript_lengths)
lengths <- as.data.frame(lengths) #convert to data frame
geneInfo$Length<-lengths$transcript_lengths[match(geneInfo$gene_id, lengths$seqnames)] #Add in length to geneInfoFormat GO terms to remove dashes and quotes and separate by semicolons (replace , with ;) in Annotation.GO.ID column
geneInfo$Annotation.GO.ID <- gsub(",", ";", geneInfo$Annotation.GO.ID)
geneInfo$Annotation.GO.ID <- gsub('"', "", geneInfo$Annotation.GO.ID)
geneInfo$Annotation.GO.ID <- gsub("-", NA, geneInfo$Annotation.GO.ID)#Calcification - upregulation (GO terms associated with WGCNA modules demonstrating positive expression in net calcification
#Run the GOSeq function by module color to test for GO term enrichment. Due to high number of enriched GO terms, I am using an adjusted p-value threshold of <0.0001 and using `rrvgo` package to reduce redundancy in list of GO terms.
##Calcification
### Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)### Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)### Generate vector with names in just the module we are analyzing
ID.vector <- geneInfo%>%
filter(moduleColor==c("blue","pink","magenta","lightcyan","purple","brown"))%>%
#get_rows(.data[[module]]))%>%
pull(gene_id)## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `==...`.
## Caused by warning in `moduleColor == c("blue", "pink", "magenta", "lightcyan", "purple", "brown")`:
## ! longer object length is not a multiple of shorter object length##Get a list of GO Terms for whole dataset
GO.terms_all <- geneInfo%>%
#filter(moduleColor==c("blue","pink","magenta","lightcyan","purple","brown"))%>%
#filter(get_rows(.data[[module]]))%>%
dplyr::select(gene_id, Annotation.GO.ID)
##Get a list of GO Terms for each module
GO.terms <- geneInfo%>%
filter(moduleColor==c("blue","pink","magenta","lightcyan","purple","brown"))%>%
#filter(get_rows(.data[[module]]))%>%
dplyr::select(gene_id, Annotation.GO.ID)## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `==...`.
## Caused by warning in `moduleColor == c("blue", "pink", "magenta", "lightcyan", "purple", "brown")`:
## ! longer object length is not a multiple of shorter object lengthhead(GO.terms)## gene_id
## 1 Pocillopora_acuta_HIv2___RNAseq.g13824.t1
## 2 Pocillopora_acuta_HIv2___RNAseq.g5347.t1
## 3 Pocillopora_acuta_HIv2___RNAseq.g16715.t1
## 4 Pocillopora_acuta_HIv2___RNAseq.g3551.t1
## 5 Pocillopora_acuta_HIv2___RNAseq.g8011.t1
## 6 Pocillopora_acuta_HIv2___RNAseq.g2726.t1
## Annotation.GO.ID
## 1 GO:0003674;GO:0003824;GO:0004064;GO:0004089;GO:0005575;GO:0005622;GO:0005623;GO:0005737;GO:0005829;GO:0006810;GO:0006811;GO:0006820;GO:0007154;GO:0007165;GO:0007166;GO:0008150;GO:0009987;GO:0010033;GO:0015701;GO:0015711;GO:0016787;GO:0016788;GO:0016829;GO:0016835;GO:0016836;GO:0019221;GO:0023052;GO:0034097;GO:0035722;GO:0042221;GO:0044424;GO:0044444;GO:0044464;GO:0050789;GO:0050794;GO:0050896;GO:0051179;GO:0051234;GO:0051716;GO:0052689;GO:0065007;GO:0070671;GO:0070887;GO:0071310;GO:0071345;GO:0071349;GO:0071702
## 2 GO:0000323;GO:0002376;GO:0002474;GO:0002478;GO:0002495;GO:0002504;GO:0003674;GO:0003824;GO:0005575;GO:0005576;GO:0005622;GO:0005623;GO:0005737;GO:0005764;GO:0005773;GO:0005775;GO:0005829;GO:0006950;GO:0006952;GO:0006955;GO:0007154;GO:0007165;GO:0007166;GO:0008150;GO:0008152;GO:0008285;GO:0009987;GO:0010033;GO:0015036;GO:0016491;GO:0016651;GO:0016667;GO:0016668;GO:0019221;GO:0019882;GO:0019884;GO:0019886;GO:0023052;GO:0030054;GO:0031647;GO:0031974;GO:0034097;GO:0034341;GO:0042127;GO:0042221;GO:0042590;GO:0043202;GO:0043226;GO:0043227;GO:0043229;GO:0043231;GO:0043233;GO:0044422;GO:0044424;GO:0044437;GO:0044444;GO:0044446;GO:0044464;GO:0045087;GO:0047134;GO:0048002;GO:0048145;GO:0048147;GO:0048519;GO:0048523;GO:0050789;GO:0050794;GO:0050821;GO:0050896;GO:0051716;GO:0055114;GO:0060333;GO:0065007;GO:0065008;GO:0070013;GO:0070887;GO:0071310;GO:0071345;GO:0071346
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 GO:0000226;GO:0001578;GO:0003341;GO:0006928;GO:0006996;GO:0007010;GO:0007017;GO:0007018;GO:0007275;GO:0007368;GO:0007389;GO:0007507;GO:0008150;GO:0009799;GO:0009855;GO:0009987;GO:0016043;GO:0022607;GO:0030030;GO:0030031;GO:0032501;GO:0032502;GO:0034622;GO:0035082;GO:0043933;GO:0044085;GO:0044458;GO:0044782;GO:0048513;GO:0048731;GO:0048856;GO:0060271;GO:0065003;GO:0070286;GO:0070925;GO:0071840;GO:0072359;GO:0120031;GO:0120036##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms$Annotation.GO.ID), ";")
split2 <- data.frame(v1 = rep.int(GO.terms$gene, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "Annotation.GO.ID")
GO.terms <- split2
dim(GO.terms)## [1] 48204 2##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms_all$Annotation.GO.ID), ";")
split2 <- data.frame(v1 = rep.int(GO.terms_all$gene, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "Annotation.GO.ID")
GO.terms_all <- split2
dim(GO.terms_all)## [1] 705483 2##Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector = as.integer(ALL.vector %in% ID.vector)
names(gene.vector)<-ALL.vector#set names
#weight gene vector by bias for length of gene
pwf <- nullp(gene.vector, ID.vector, bias.data=LENGTH.vector) ## Warning in pcls(G): initial point very close to some inequality constraints
#run goseq using Wallenius method for all categories of GO terms
GO.wall <- goseq(pwf, ID.vector, gene2cat=GO.terms, test.cats=c("GO:BP", "GO:MF", "GO:CC"), method="Wallenius", use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columnsGO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"#adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method="BH") #add adjusted p-values#Filtering for p < 0.00001
GO <- GO %>%
dplyr::filter(bh_adjust<0.00001) %>%
dplyr::arrange(., ontology, bh_adjust)
#Write file of results
write.csv(GO, file = "../../../output/WGCNA/goseq_pattern_calcification.csv")module_vector <- c(levels(geneInfo$moduleColor))
ontologies <- c("BP", "MF")go_results <- read.csv("../../../output/WGCNA/goseq_pattern_calcification.csv")
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0000375 0 1 13
## 2 2 GO:0000377 0 1 13
## 3 3 GO:0000398 0 1 13
## 4 4 GO:0001505 0 1 16
## 5 5 GO:0001568 0 1 11
## 6 6 GO:0001655 0 1 14
## numInCat
## 1 13
## 2 13
## 3 13
## 4 16
## 5 11
## 6 14
## term
## 1 RNA splicing, via transesterification reactions
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 mRNA splicing, via spliceosome
## 4 regulation of neurotransmitter levels
## 5 blood vessel development
## 6 urogenital system development
## ontology bh_adjust
## 1 BP 0
## 2 BP 0
## 3 BP 0
## 4 BP 0
## 5 BP 0
## 6 BP 0go_results <- go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
filter(numInCat>10)%>%
arrange(., bh_adjust)
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0000375 0 1 13
## 2 2 GO:0000377 0 1 13
## 3 3 GO:0000398 0 1 13
## 4 4 GO:0001505 0 1 16
## 5 5 GO:0001568 0 1 11
## 6 6 GO:0001655 0 1 14
## numInCat
## 1 13
## 2 13
## 3 13
## 4 16
## 5 11
## 6 14
## term
## 1 RNA splicing, via transesterification reactions
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 mRNA splicing, via spliceosome
## 4 regulation of neurotransmitter levels
## 5 blood vessel development
## 6 urogenital system development
## ontology bh_adjust
## 1 BP 0
## 2 BP 0
## 3 BP 0
## 4 BP 0
## 5 BP 0
## 6 BP 0library(rrvgo)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")## ## preparing gene to GO mapping data...## preparing IC data...#calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
head(reducedTerms)## go cluster parent score size
## GO:0098657 GO:0098657 61 GO:0098657 Inf 0
## GO:0032259 GO:0032259 50 GO:0032259 Inf 107
## GO:0009987 GO:0009987 36 GO:0009987 Inf 0
## GO:0007568 GO:0007568 26 GO:0007568 Inf 19
## GO:0007163 GO:0007163 21 GO:0007163 Inf 27
## GO:0005975 GO:0005975 8 GO:0005975 Inf 161
## term
## GO:0098657 import into cell
## GO:0032259 methylation
## GO:0009987 cellular process
## GO:0007568 aging
## GO:0007163 establishment or maintenance of cell polarity
## GO:0005975 carbohydrate metabolic process
## parentTerm termUniqueness
## GO:0098657 import into cell 0.9673904
## GO:0032259 methylation 0.9627408
## GO:0009987 cellular process 0.9690624
## GO:0007568 aging 0.9482864
## GO:0007163 establishment or maintenance of cell polarity 0.9848672
## GO:0005975 carbohydrate metabolic process 0.9538752
## termUniquenessWithinCluster termDispensability
## GO:0098657 1 0
## GO:0032259 1 0
## GO:0009987 1 0
## GO:0007568 1 0
## GO:0007163 1 0
## GO:0005975 1 0#keep only the goterms from the reduced list
go_results <- go_results%>%
filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm <- reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]
go_results <- go_results %>% mutate(Factor = "Up")
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0000375 0 1 13
## 2 2 GO:0000377 0 1 13
## 3 3 GO:0000398 0 1 13
## 4 4 GO:0001505 0 1 16
## 5 5 GO:0001568 0 1 11
## 6 9 GO:0001775 0 1 16
## numInCat
## 1 13
## 2 13
## 3 13
## 4 16
## 5 11
## 6 16
## term
## 1 RNA splicing, via transesterification reactions
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 mRNA splicing, via spliceosome
## 4 regulation of neurotransmitter levels
## 5 blood vessel development
## 6 cell activation
## ontology bh_adjust
## 1 BP 0
## 2 BP 0
## 3 BP 0
## 4 BP 0
## 5 BP 0
## 6 BP 0
## ParentTerm
## 1 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 4 regulation of hormone secretion
## 5 gliogenesis
## 6 leukocyte activation
## Factor
## 1 Up
## 2 Up
## 3 Up
## 4 Up
## 5 Up
## 6 Upwrite.csv(go_results, "../../../output/WGCNA/goseq_pattern_calcification_filtered.csv")#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot_calcification <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot_calcification
#Calcification - downregulation (GO terms associated with WGCNA modules demonstrating negative expression in net calcification
### Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)### Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)### Generate vector with names in just the module we are analyzing
ID.vector <- geneInfo%>%
filter(moduleColor==c("black","red"))%>%
#get_rows(.data[[module]]))%>%
pull(gene_id)## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `moduleColor == c("black", "red")`.
## Caused by warning in `moduleColor == c("black", "red")`:
## ! longer object length is not a multiple of shorter object length##Get a list of GO Terms for each module
GO.terms <- geneInfo%>%
filter(moduleColor==c("black","red"))%>%
#filter(get_rows(.data[[module]]))%>%
dplyr::select(gene_id, Annotation.GO.ID)## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `moduleColor == c("black", "red")`.
## Caused by warning in `moduleColor == c("black", "red")`:
## ! longer object length is not a multiple of shorter object lengthhead(GO.terms)## gene_id
## 1 Pocillopora_acuta_HIv2___TS.g11738.t1
## 2 Pocillopora_acuta_HIv2___RNAseq.g11877.t2b
## 3 Pocillopora_acuta_HIv2___RNAseq.g11040.t1
## 4 Pocillopora_acuta_HIv2___RNAseq.g11055.t1
## 5 Pocillopora_acuta_HIv2___RNAseq.g7761.t1
## 6 Pocillopora_acuta_HIv2___RNAseq.g7483.t2c
## Annotation.GO.ID
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 GO:0003674;GO:0003824;GO:0005575;GO:0005622;GO:0005623;GO:0005737;GO:0005783;GO:0006464;GO:0006497;GO:0006807;GO:0008150;GO:0008152;GO:0009058;GO:0009059;GO:0009987;GO:0012505;GO:0016409;GO:0016417;GO:0016740;GO:0016746;GO:0016747;GO:0018345;GO:0019538;GO:0019706;GO:0019707;GO:0034645;GO:0036211;GO:0042157;GO:0042158;GO:0043170;GO:0043226;GO:0043227;GO:0043229;GO:0043231;GO:0043412;GO:0043543;GO:0044237;GO:0044238;GO:0044249;GO:0044260;GO:0044267;GO:0044424;GO:0044444;GO:0044464;GO:0071704;GO:0140096;GO:1901564;GO:1901566;GO:1901576
## 6 <NA>##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms$Annotation.GO.ID), ";")
split2 <- data.frame(v1 = rep.int(GO.terms$gene, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "Annotation.GO.ID")
GO.terms <- split2
dim(GO.terms)## [1] 28453 2##Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector=as.integer(ALL.vector %in% ID.vector)
names(gene.vector) <- ALL.vector#set names
#weight gene vector by bias for length of gene
pwf <- nullp(gene.vector, ID.vector, bias.data=LENGTH.vector) 
#run goseq using Wallenius method for all categories of GO terms
GO.wall<-goseq(pwf, ID.vector, gene2cat=GO.terms, test.cats=c("GO:BP", "GO:MF", "GO:CC"), method="Wallenius", use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columnsGO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"#adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method="BH") #add adjusted p-values#Filtering for p < 0.01
GO <- GO %>%
dplyr::filter(bh_adjust<0.00001) %>%
dplyr::arrange(., ontology, bh_adjust)
#Write file of results
write.csv(GO, file = "../../../output/WGCNA/goseq_pattern_calcification_down.csv")module_vector <- c(levels(geneInfo$moduleColor))
ontologies <- c("BP", "MF")go_results <- read.csv("../../../output/WGCNA/goseq_pattern_calcification_down.csv")go_results<-go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
filter(numInCat>10)%>%
arrange(., bh_adjust)
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 2 GO:0000226 0 1 11
## 2 3 GO:0000278 0 1 21
## 3 4 GO:0000280 0 1 11
## 4 5 GO:0000902 0 1 16
## 5 6 GO:0000904 0 1 12
## 6 7 GO:0001101 0 1 13
## numInCat term ontology bh_adjust
## 1 11 microtubule cytoskeleton organization BP 0
## 2 21 mitotic cell cycle BP 0
## 3 11 nuclear division BP 0
## 4 16 cell morphogenesis BP 0
## 5 12 cell morphogenesis involved in differentiation BP 0
## 6 13 response to acid chemical BP 0library(rrvgo)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")## preparing gene to GO mapping data...## preparing IC data... #calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
head(reducedTerms)## go cluster parent score size
## GO:0051301 GO:0051301 43 GO:0051301 Inf 162
## GO:0043603 GO:0043603 41 GO:0043603 Inf 0
## GO:0030097 GO:0030097 37 GO:0030097 Inf 0
## GO:0007610 GO:0007610 13 GO:0007610 Inf 18
## GO:0090066 GO:0090066 45 GO:0090066 Inf 0
## GO:0001101 GO:0001101 5 GO:0001101 Inf 0
## term
## GO:0051301 cell division
## GO:0043603 amide metabolic process
## GO:0030097 hemopoiesis
## GO:0007610 behavior
## GO:0090066 regulation of anatomical structure size
## GO:0001101 response to acid chemical
## parentTerm termUniqueness
## GO:0051301 cell division 0.9834650
## GO:0043603 amide metabolic process 0.9348500
## GO:0030097 hemopoiesis 0.9314200
## GO:0007610 behavior 0.9518225
## GO:0090066 regulation of anatomical structure size 0.9279425
## GO:0001101 response to acid chemical 0.9681175
## termUniquenessWithinCluster termDispensability
## GO:0051301 1.0000000 0
## GO:0043603 1.0000000 0
## GO:0030097 1.0000000 0
## GO:0007610 0.6508000 0
## GO:0090066 0.6095000 0
## GO:0001101 0.6048667 0#keep only the goterms from the reduced list
go_results<-go_results%>%
filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm<-reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]
go_results <- go_results %>% mutate(Factor = "Down")
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 2 GO:0000226 0 1 11
## 2 3 GO:0000278 0 1 21
## 3 4 GO:0000280 0 1 11
## 4 5 GO:0000902 0 1 16
## 5 6 GO:0000904 0 1 12
## 6 7 GO:0001101 0 1 13
## numInCat term ontology bh_adjust
## 1 11 microtubule cytoskeleton organization BP 0
## 2 21 mitotic cell cycle BP 0
## 3 11 nuclear division BP 0
## 4 16 cell morphogenesis BP 0
## 5 12 cell morphogenesis involved in differentiation BP 0
## 6 13 response to acid chemical BP 0
## ParentTerm Factor
## 1 microtubule-based process Down
## 2 negative regulation of cell cycle Down
## 3 endomembrane system organization Down
## 4 regionalization Down
## 5 regionalization Down
## 6 response to acid chemical Downwrite.csv(go_results, "../../../output/WGCNA/goseq_pattern_calcification_down_filtered.csv")#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot_calcification_down <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot_calcification_down
#Frequency of GO terms
cal_up_terms <-read.csv("../../../output/WGCNA/goseq_pattern_calcification_filtered.csv")
cal_down_terms <-read.csv("../../../output/WGCNA/goseq_pattern_calcification_down_filtered.csv")cal_biomin_terms <-read.csv("../../../output/Biomineralization_goterms.csv")
head(cal_biomin_terms)## Factor X.1 X GOterm over_represented_pvalue under_represented_pvalue
## 1 Biomin 110 158 GO:0030048 1 0.9417955
## 2 Biomin 135 192 GO:0070252 1 0.9417955
## 3 Biomin 397 563 GO:0006928 1 0.8194394
## 4 Biomin 449 627 GO:0030334 1 0.9494130
## 5 Biomin 450 628 GO:0030335 1 0.9494130
## 6 Biomin 451 629 GO:0030336 1 0.9494130
## numDEInCat numInCat term ontology
## 1 0 1 actin filament-based movement BP
## 2 0 1 actin-mediated cell contraction BP
## 3 0 4 movement of cell or subcellular component BP
## 4 0 1 regulation of cell migration BP
## 5 0 1 positive regulation of cell migration BP
## 6 0 1 negative regulation of cell migration BP
## bh_adjust ParentTerm
## 1 1 actin-mediated cell contraction
## 2 1 actin-mediated cell contraction
## 3 1 actin-mediated cell contraction
## 4 1 actin-mediated cell contraction
## 5 1 actin-mediated cell contraction
## 6 1 actin-mediated cell contractionhead(cal_up_terms)## X.1 X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 1 GO:0000375 0 1 13
## 2 2 2 GO:0000377 0 1 13
## 3 3 3 GO:0000398 0 1 13
## 4 4 4 GO:0001505 0 1 16
## 5 5 5 GO:0001568 0 1 11
## 6 6 9 GO:0001775 0 1 16
## numInCat
## 1 13
## 2 13
## 3 13
## 4 16
## 5 11
## 6 16
## term
## 1 RNA splicing, via transesterification reactions
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 mRNA splicing, via spliceosome
## 4 regulation of neurotransmitter levels
## 5 blood vessel development
## 6 cell activation
## ontology bh_adjust
## 1 BP 0
## 2 BP 0
## 3 BP 0
## 4 BP 0
## 5 BP 0
## 6 BP 0
## ParentTerm
## 1 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 2 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 3 RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
## 4 regulation of hormone secretion
## 5 gliogenesis
## 6 leukocyte activation
## Factor
## 1 Up
## 2 Up
## 3 Up
## 4 Up
## 5 Up
## 6 Uphead(cal_down_terms)## X.1 X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 2 GO:0000226 0 1 11
## 2 2 3 GO:0000278 0 1 21
## 3 3 4 GO:0000280 0 1 11
## 4 4 5 GO:0000902 0 1 16
## 5 5 6 GO:0000904 0 1 12
## 6 6 7 GO:0001101 0 1 13
## numInCat term ontology bh_adjust
## 1 11 microtubule cytoskeleton organization BP 0
## 2 21 mitotic cell cycle BP 0
## 3 11 nuclear division BP 0
## 4 16 cell morphogenesis BP 0
## 5 12 cell morphogenesis involved in differentiation BP 0
## 6 13 response to acid chemical BP 0
## ParentTerm Factor
## 1 microtubule-based process Down
## 2 negative regulation of cell cycle Down
## 3 endomembrane system organization Down
## 4 regionalization Down
## 5 regionalization Down
## 6 response to acid chemical Downall_terms<- merge(cal_up_terms,cal_down_terms, by=c("Factor","GOterm","X.1","X","GOterm","over_represented_pvalue","under_represented_pvalue","numDEInCat","numInCat","term","ontology","bh_adjust","ParentTerm"),all=T)
all_terms<- merge(all_terms,cal_biomin_terms, by=c("Factor","GOterm","X.1","X","GOterm","over_represented_pvalue","under_represented_pvalue","numDEInCat","numInCat","term","ontology","bh_adjust","ParentTerm"),all=T)
all_terms$GOterm<-as.factor(all_terms$GOterm)
head(all_terms)## Factor GOterm X.1 X over_represented_pvalue under_represented_pvalue
## 1 Biomin GO:0000003 748 1054 1 0.7918886
## 2 Biomin GO:0000041 335 479 1 0.9545025
## 3 Biomin GO:0000122 295 438 1 0.9689632
## 4 Biomin GO:0000132 97 143 1 0.9417955
## 5 Biomin GO:0000226 390 556 1 0.8588726
## 6 Biomin GO:0000278 595 812 1 0.9325796
## numDEInCat numInCat term
## 1 0 5 reproduction
## 2 0 1 transition metal ion transport
## 3 0 1 negative regulation of transcription by RNA polymerase II
## 4 0 1 establishment of mitotic spindle orientation
## 5 0 3 microtubule cytoskeleton organization
## 6 0 2 mitotic cell cycle
## ontology bh_adjust
## 1 BP 1
## 2 BP 1
## 3 BP 1
## 4 BP 1
## 5 BP 1
## 6 BP 1
## ParentTerm
## 1 gonad development
## 2 divalent metal ion transport
## 3 negative regulation of nucleobase-containing compound metabolic process
## 4 establishment of mitotic spindle orientation
## 5 microtubule-based process
## 6 microtubule cytoskeleton organization involved in mitosisall_terms_merged <- all_terms %>%
group_by(GOterm) %>%
dplyr::summarise(
ParentTerm = paste(unique(ParentTerm), collapse = ", "),
Factor = paste(unique(Factor), collapse = ", ")
)
head(all_terms_merged)## # A tibble: 6 × 3
## GOterm ParentTerm Factor
## <fct> <chr> <chr>
## 1 GO:0000003 gonad development, development of primary sexual characteri… Biomi…
## 2 GO:0000041 divalent metal ion transport Biomin
## 3 GO:0000122 negative regulation of nucleobase-containing compound metab… Biomi…
## 4 GO:0000132 establishment of mitotic spindle orientation Biomin
## 5 GO:0000226 microtubule-based process, microtubule-based movement Biomi…
## 6 GO:0000278 microtubule cytoskeleton organization involved in mitosis, … Biomi…write.csv(all_terms_merged, "../../../output/WGCNA/Merged_GOterms_Factor_ParentTerm.csv")cal_freq_terms <-read.csv("../../../output/WGCNA/goseq_pattern_calcification_all.csv")
head(cal_freq_terms)## ParentTerm Calcification.direction
## 1 negative regulation of organelle organization Up
## 2 negative regulation of protein modification process Up
## 3 anion transport Up
## 4 sensory organ morphogenesis Up
## 5 cytokine-mediated signaling pathway Up
## 6 biological regulation Up
## Number.of.terms Frequency.of.shared.biomin.GOterms
## 1 39 12
## 2 26 11
## 3 23 6
## 4 21 3
## 5 20 2
## 6 19 9
## Proportion.of.shared.GO.terms.with.biomineralization.genes
## 1 0.3076923
## 2 0.4230769
## 3 0.2608696
## 4 0.1428571
## 5 0.1000000
## 6 0.4736842
## Percentage.of.shared.GO.terms.with.biomineralization.genes
## 1 30.76923
## 2 42.30769
## 3 26.08696
## 4 14.28571
## 5 10.00000
## 6 47.36842##Frequency >10 GO terms upregulation
cal_freq_terms_filtered_up <- cal_freq_terms %>%
filter(Number.of.terms>=10) %>%
filter(Calcification.direction=="Up")###Figure
#counts$Direction.of.flat.origin<- factor(counts$Direction.of.flat.origin, levels =c("up","no pattern","down"))
#counts$Module<- factor(counts$Module, levels=c("Blue","Brown","Greenyellow","Cyan","Pink","Magenta","Lightcyan","Midnight blue","Purple","Turquiose","Red","Black"))
freq_fig_up<-ggplot(cal_freq_terms_filtered_up, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes,group=1))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,40))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(0, 100))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95,size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.freq_fig_up
##Frequency >10 GO terms downregulation
cal_freq_terms_filtered_down <- cal_freq_terms %>%
filter(Number.of.terms>=10) %>%
filter(Calcification.direction=="Down")###Figure
freq_fig_down<-ggplot(cal_freq_terms_filtered_down, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
#geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,40))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#d1e5f0","#92c5de","#4393c3","#2166ac"), na.value = "grey98",limits = c(0, 100))+
#scale_fill_gradientn(colours=c("#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(10, 40))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95, size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_down
compare_figs<-cowplot::plot_grid(freq_fig_up, freq_fig_down, nrow=2, align="v")
compare_figs
##Frequency all GO terms upregulation
cal_freq_terms_filtered_up_all <- cal_freq_terms %>%
filter(Calcification.direction=="Up")
head(cal_freq_terms_filtered_up_all)## ParentTerm Calcification.direction
## 1 negative regulation of organelle organization Up
## 2 negative regulation of protein modification process Up
## 3 anion transport Up
## 4 sensory organ morphogenesis Up
## 5 cytokine-mediated signaling pathway Up
## 6 biological regulation Up
## Number.of.terms Frequency.of.shared.biomin.GOterms
## 1 39 12
## 2 26 11
## 3 23 6
## 4 21 3
## 5 20 2
## 6 19 9
## Proportion.of.shared.GO.terms.with.biomineralization.genes
## 1 0.3076923
## 2 0.4230769
## 3 0.2608696
## 4 0.1428571
## 5 0.1000000
## 6 0.4736842
## Percentage.of.shared.GO.terms.with.biomineralization.genes
## 1 30.76923
## 2 42.30769
## 3 26.08696
## 4 14.28571
## 5 10.00000
## 6 47.36842###Figure
#counts$Direction.of.flat.origin<- factor(counts$Direction.of.flat.origin, levels =c("up","no pattern","down"))
#counts$Module<- factor(counts$Module, levels=c("Blue","Brown","Greenyellow","Cyan","Pink","Magenta","Lightcyan","Midnight blue","Purple","Turquiose","Red","Black"))
freq_fig_up<-ggplot(cal_freq_terms_filtered_up_all, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes,group=1))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,40))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(0, 100))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95,size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_up
##Frequency all GO terms downregulation
cal_freq_terms_filtered_down_all <- cal_freq_terms %>%
filter(Calcification.direction=="Down")
cal_freq_terms_filtered_down_all## ParentTerm
## 1 aging
## 2 axon guidance
## 3 biological regulation
## 4 carbohydrate derivative biosynthetic process
## 5 carbohydrate metabolic process
## 6 cation transport
## 7 cellular lipid metabolic process
## 8 cellular process
## 9 cellular response to abiotic stimulus
## 10 cellular response to external stimulus
## 11 cellular response to hormone stimulus
## 12 detection of stimulus
## 13 development of primary sexual characteristics
## 14 developmental growth
## 15 DNA repair
## 16 drug metabolic process
## 17 embryonic organ development
## 18 establishment or maintenance of cell polarity
## 19 exocytosis
## 20 gastrulation
## 21 head development
## 22 histone modification
## 23 import into cell
## 24 leukocyte activation
## 25 locomotory behavior
## 26 macromolecule metabolic process
## 27 memory
## 28 metabolic process
## 29 morphogenesis of embryonic epithelium
## 30 negative regulation of catabolic process
## 31 negative regulation of cell development
## 32 negative regulation of mitotic cell cycle
## 33 negative regulation of neurogenesis
## 34 negative regulation of transcription by RNA polymerase II
## 35 nucleotide metabolic process
## 36 organonitrogen compound metabolic process
## 37 oxidation-reduction process
## 38 positive regulation of cell motility
## 39 positive regulation of defense response
## 40 positive regulation of immune response
## 41 positive regulation of kinase activity
## 42 positive regulation of multi-organism process
## 43 post-embryonic animal morphogenesis
## 44 protein localization to organelle
## 45 protein modification by small protein conjugation or removal
## 46 purine ribonucleotide metabolic process
## 47 regulation of actin cytoskeleton organization
## 48 regulation of cell population proliferation
## 49 regulation of cell size
## 50 regulation of cell-cell adhesion
## 51 regulation of gene expression, epigenetic
## 52 regulation of ion transport
## 53 regulation of microtubule-based process
## 54 regulation of neuron death
## 55 regulation of peptide secretion
## 56 regulation of Ras protein signal transduction
## 57 response to ionizing radiation
## 58 response to xenobiotic stimulus
## 59 ribose phosphate metabolic process
## 60 rRNA metabolic process
## 61 transcription by RNA polymerase II
## 62 transmembrane receptor protein tyrosine kinase signaling pathway
## 63 tube development
## Calcification.direction Number.of.terms Frequency.of.shared.biomin.GOterms
## 1 Down 12 11
## 2 Down 7 7
## 3 Down 16 9
## 4 Down 2 1
## 5 Down 1 0
## 6 Down 14 14
## 7 Down 3 1
## 8 Down 1 1
## 9 Down 2 0
## 10 Down 5 5
## 11 Down 4 4
## 12 Down 6 3
## 13 Down 10 8
## 14 Down 8 6
## 15 Down 8 5
## 16 Down 1 1
## 17 Down 19 19
## 18 Down 1 1
## 19 Down 9 1
## 20 Down 7 7
## 21 Down 2 2
## 22 Down 39 27
## 23 Down 1 1
## 24 Down 2 0
## 25 Down 2 2
## 26 Down 3 3
## 27 Down 9 5
## 28 Down 6 6
## 29 Down 7 7
## 30 Down 13 12
## 31 Down 10 10
## 32 Down 10 7
## 33 Down 13 11
## 34 Down 17 17
## 35 Down 7 7
## 36 Down 6 4
## 37 Down 1 1
## 38 Down 8 7
## 39 Down 7 4
## 40 Down 9 1
## 41 Down 20 17
## 42 Down 10 6
## 43 Down 19 19
## 44 Down 14 9
## 45 Down 6 4
## 46 Down 13 7
## 47 Down 4 2
## 48 Down 4 1
## 49 Down 20 18
## 50 Down 2 0
## 51 Down 27 26
## 52 Down 17 7
## 53 Down 3 2
## 54 Down 13 5
## 55 Down 8 0
## 56 Down 7 6
## 57 Down 7 5
## 58 Down 25 19
## 59 Down 22 12
## 60 Down 6 1
## 61 Down 14 1
## 62 Down 23 21
## 63 Down 3 3
## Proportion.of.shared.GO.terms.with.biomineralization.genes
## 1 0.91666667
## 2 1.00000000
## 3 0.56250000
## 4 0.50000000
## 5 0.00000000
## 6 1.00000000
## 7 0.33333333
## 8 1.00000000
## 9 0.00000000
## 10 1.00000000
## 11 1.00000000
## 12 0.50000000
## 13 0.80000000
## 14 0.75000000
## 15 0.62500000
## 16 1.00000000
## 17 1.00000000
## 18 1.00000000
## 19 0.11111111
## 20 1.00000000
## 21 1.00000000
## 22 0.69230769
## 23 1.00000000
## 24 0.00000000
## 25 1.00000000
## 26 1.00000000
## 27 0.55555556
## 28 1.00000000
## 29 1.00000000
## 30 0.92307692
## 31 1.00000000
## 32 0.70000000
## 33 0.84615385
## 34 1.00000000
## 35 1.00000000
## 36 0.66666667
## 37 1.00000000
## 38 0.87500000
## 39 0.57142857
## 40 0.11111111
## 41 0.85000000
## 42 0.60000000
## 43 1.00000000
## 44 0.64285714
## 45 0.66666667
## 46 0.53846154
## 47 0.50000000
## 48 0.25000000
## 49 0.90000000
## 50 0.00000000
## 51 0.96296296
## 52 0.41176471
## 53 0.66666667
## 54 0.38461538
## 55 0.00000000
## 56 0.85714286
## 57 0.71428571
## 58 0.76000000
## 59 0.54545455
## 60 0.16666667
## 61 0.07142857
## 62 0.91304348
## 63 1.00000000
## Percentage.of.shared.GO.terms.with.biomineralization.genes
## 1 91.666667
## 2 100.000000
## 3 56.250000
## 4 50.000000
## 5 0.000000
## 6 100.000000
## 7 33.333333
## 8 100.000000
## 9 0.000000
## 10 100.000000
## 11 100.000000
## 12 50.000000
## 13 80.000000
## 14 75.000000
## 15 62.500000
## 16 100.000000
## 17 100.000000
## 18 100.000000
## 19 11.111111
## 20 100.000000
## 21 100.000000
## 22 69.230769
## 23 100.000000
## 24 0.000000
## 25 100.000000
## 26 100.000000
## 27 55.555556
## 28 100.000000
## 29 100.000000
## 30 92.307692
## 31 100.000000
## 32 70.000000
## 33 84.615385
## 34 100.000000
## 35 100.000000
## 36 66.666667
## 37 100.000000
## 38 87.500000
## 39 57.142857
## 40 11.111111
## 41 85.000000
## 42 60.000000
## 43 100.000000
## 44 64.285714
## 45 66.666667
## 46 53.846154
## 47 50.000000
## 48 25.000000
## 49 90.000000
## 50 0.000000
## 51 96.296296
## 52 41.176471
## 53 66.666667
## 54 38.461538
## 55 0.000000
## 56 85.714286
## 57 71.428571
## 58 76.000000
## 59 54.545455
## 60 16.666667
## 61 7.142857
## 62 91.304348
## 63 100.000000###Figure
freq_fig_down<-ggplot(cal_freq_terms_filtered_down_all, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
#geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,40))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#d1e5f0","#92c5de","#4393c3","#2166ac"), na.value = "grey98",limits = c(0, 100))+
#scale_fill_gradientn(colours=c("#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(10, 40))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95, size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_down
compare_figs_all<-cowplot::plot_grid(freq_fig_up, freq_fig_down, ncol=2, align="h")
compare_figs_all
biomin_mod <-read.csv("../../../output/WGCNA/Biomineralization_toolkit_by_module.csv")
head(biomin_mod)## X Pocillopora_acuta_best_hit accessionnumber.geneID
## 1 225 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 XP_022804785.1
## 2 608 Pocillopora_acuta_HIv2___RNAseq.g11609.t1 P33_g8985
## 3 1253 Pocillopora_acuta_HIv2___RNAseq.g13824.t1 Gene:g27814
## 4 1427 Pocillopora_acuta_HIv2___RNAseq.g14505.t1 Gene:g9094
## 5 1466 Pocillopora_acuta_HIv2___RNAseq.g14663.t1a PFX27832.1
## 6 1642 Pocillopora_acuta_HIv2___RNAseq.g15280.t1 AJQ31790.1
## definition
## 1 thioredoxin reductase 1, cytoplasmic-like [Stylophora pistillata]
## 2 Flagellar associated protein
## 3 Annotated: carbonic anhydrase (STPCA2-2)
## 4 Annotated: Actin
## 5 Poly [ADP-ribose] polymerase 11 [Stylophora pistillata]
## 6 solute carrier family 4 member gamma [Stylophora pistillata]
## Ref substanceBXH
## 1 Peled et al., 2020 Pocillopora_acuta_HIv2___RNAseq.g10093.t2
## 2 Drake et al., 2013 Pocillopora_acuta_HIv2___RNAseq.g11609.t1
## 3 Mummadisetti et al., 2021 Pocillopora_acuta_HIv2___RNAseq.g13824.t1
## 4 Mummadisetti et al., 2021 Pocillopora_acuta_HIv2___RNAseq.g14505.t1
## 5 Peled et al., 2020 Pocillopora_acuta_HIv2___RNAseq.g14663.t1a
## 6 Zoccola et al., 2015 Pocillopora_acuta_HIv2___RNAseq.g15280.t1
## geneSymbol moduleColor
## 1 Pocillopora_acuta_HIv2___Sc0000021 blue
## 2 Pocillopora_acuta_HIv2___Sc0000013 turquoise
## 3 Pocillopora_acuta_HIv2___Sc0000005 blue
## 4 Pocillopora_acuta_HIv2___xfSc0000036 turquoise
## 5 Pocillopora_acuta_HIv2___Sc0000020 turquoise
## 6 Pocillopora_acuta_HIv2___Sc0000006 pink
## GO.terms
## 1 GO:0000003,GO:0000302,GO:0000305,GO:0001650,GO:0001704,GO:0001707,GO:0001887,GO:0001890,GO:0003006,GO:0003674,GO:0003824,GO:0004791,GO:0005488,GO:0005515,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005654,GO:0005730,GO:0005737,GO:0005739,GO:0005783,GO:0005829,GO:0006082,GO:0006139,GO:0006518,GO:0006520,GO:0006575,GO:0006725,GO:0006732,GO:0006733,GO:0006739,GO:0006749,GO:0006753,GO:0006790,GO:0006793,GO:0006796,GO:0006807,GO:0006950,GO:0006979,GO:0007154,GO:0007165,GO:0007275,GO:0007369,GO:0007498,GO:0008150,GO:0008152,GO:0008283,GO:0009056,GO:0009069,GO:0009117,GO:0009611,GO:0009628,GO:0009636,GO:0009653,GO:0009790,GO:0009888,GO:0009987,GO:0010035,GO:0010038,GO:0010269,GO:0010941,GO:0010942,GO:0012505,GO:0015036,GO:0015949,GO:0016043,GO:0016174,GO:0016209,GO:0016259,GO:0016491,GO:0016651,GO:0016667,GO:0016668,GO:0016999,GO:0017001,GO:0017144,GO:0018996,GO:0019216,GO:0019222,GO:0019362,GO:0019637,GO:0019725,GO:0019752,GO:0022404,GO:0022414,GO:0022607,GO:0023052,GO:0031974,GO:0031981,GO:0032501,GO:0032502,GO:0033554,GO:0033797,GO:0034599,GO:0034641,GO:0036295,GO:0036296,GO:0036477,GO:0042221,GO:0042303,GO:0042395,GO:0042493,GO:0042537,GO:0042592,GO:0042737,GO:0042743,GO:0042744,GO:0042802,GO:0042803,GO:0043025,GO:0043167,GO:0043169,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043233,GO:0043436,GO:0043603,GO:0043933,GO:0044085,GO:0044237,GO:0044238,GO:0044248,GO:0044281,GO:0044297,GO:0044422,GO:0044424,GO:0044428,GO:0044444,GO:0044446,GO:0044452,GO:0044464,GO:0045340,GO:0045454,GO:0046483,GO:0046496,GO:0046688,GO:0046872,GO:0046914,GO:0046983,GO:0048332,GO:0048513,GO:0048518,GO:0048522,GO:0048598,GO:0048608,GO:0048646,GO:0048678,GO:0048729,GO:0048731,GO:0048856,GO:0050664,GO:0050789,GO:0050794,GO:0050896,GO:0051186,GO:0051187,GO:0051259,GO:0051262,GO:0051716,GO:0055086,GO:0055093,GO:0055114,GO:0061458,GO:0065003,GO:0065007,GO:0065008,GO:0070013,GO:0070276,GO:0070482,GO:0070887,GO:0070995,GO:0071241,GO:0071248,GO:0071280,GO:0071453,GO:0071455,GO:0071704,GO:0071840,GO:0072524,GO:0072593,GO:0080090,GO:0097237,GO:0097458,GO:0098623,GO:0098625,GO:0098626,GO:0098754,GO:0098869,GO:1901360,GO:1901564,GO:1901605,GO:1901700,GO:1990748
## 2 -
## 3 GO:0003674,GO:0003824,GO:0004064,GO:0004089,GO:0005575,GO:0005622,GO:0005623,GO:0005737,GO:0005829,GO:0006810,GO:0006811,GO:0006820,GO:0007154,GO:0007165,GO:0007166,GO:0008150,GO:0009987,GO:0010033,GO:0015701,GO:0015711,GO:0016787,GO:0016788,GO:0016829,GO:0016835,GO:0016836,GO:0019221,GO:0023052,GO:0034097,GO:0035722,GO:0042221,GO:0044424,GO:0044444,GO:0044464,GO:0050789,GO:0050794,GO:0050896,GO:0051179,GO:0051234,GO:0051716,GO:0052689,GO:0065007,GO:0070671,GO:0070887,GO:0071310,GO:0071345,GO:0071349,GO:0071702
## 4 GO:0000123,GO:0000278,GO:0000281,GO:0000910,GO:0003674,GO:0005198,GO:0005200,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005654,GO:0005737,GO:0005856,GO:0005884,GO:0005886,GO:0005912,GO:0005924,GO:0005925,GO:0006325,GO:0006338,GO:0006464,GO:0006473,GO:0006475,GO:0006807,GO:0006996,GO:0007010,GO:0007015,GO:0007049,GO:0007275,GO:0007507,GO:0007517,GO:0007519,GO:0007528,GO:0008150,GO:0008152,GO:0009653,GO:0009888,GO:0009987,GO:0010927,GO:0014706,GO:0014866,GO:0014902,GO:0014904,GO:0015629,GO:0016020,GO:0016043,GO:0016569,GO:0016570,GO:0016573,GO:0018193,GO:0018205,GO:0018393,GO:0018394,GO:0019538,GO:0022402,GO:0022607,GO:0030029,GO:0030036,GO:0030054,GO:0030055,GO:0030154,GO:0030239,GO:0031032,GO:0031248,GO:0031974,GO:0031981,GO:0032501,GO:0032502,GO:0032989,GO:0032991,GO:0034728,GO:0035267,GO:0036211,GO:0042692,GO:0043044,GO:0043170,GO:0043189,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043233,GO:0043412,GO:0043486,GO:0043543,GO:0043933,GO:0044085,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044422,GO:0044424,GO:0044428,GO:0044430,GO:0044444,GO:0044446,GO:0044451,GO:0044464,GO:0048468,GO:0048513,GO:0048646,GO:0048731,GO:0048741,GO:0048747,GO:0048856,GO:0048869,GO:0050789,GO:0050794,GO:0050803,GO:0050807,GO:0050808,GO:0051128,GO:0051146,GO:0051276,GO:0051301,GO:0055001,GO:0055002,GO:0060537,GO:0060538,GO:0061061,GO:0061640,GO:0065007,GO:0065008,GO:0070013,GO:0070161,GO:0070925,GO:0071689,GO:0071704,GO:0071824,GO:0071840,GO:0071944,GO:0072359,GO:0097433,GO:0097435,GO:0098529,GO:0099080,GO:0099081,GO:0099512,GO:0099513,GO:1901564,GO:1902493,GO:1902494,GO:1902562,GO:1903047,GO:1990234
## 5 GO:0000003,GO:0001067,GO:0001501,GO:0001568,GO:0001570,GO:0001655,GO:0001822,GO:0001944,GO:0002376,GO:0002520,GO:0003006,GO:0003674,GO:0003676,GO:0003677,GO:0003824,GO:0003950,GO:0005488,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0006464,GO:0006471,GO:0006629,GO:0006807,GO:0007154,GO:0007165,GO:0007166,GO:0007167,GO:0007169,GO:0007275,GO:0007548,GO:0008150,GO:0008152,GO:0008202,GO:0008209,GO:0008210,GO:0008406,GO:0008585,GO:0009653,GO:0009791,GO:0009887,GO:0009888,GO:0009892,GO:0009893,GO:0009894,GO:0009896,GO:0009987,GO:0010033,GO:0010171,GO:0010468,GO:0010604,GO:0010605,GO:0010629,GO:0010817,GO:0014070,GO:0016740,GO:0016757,GO:0016763,GO:0019222,GO:0019538,GO:0022414,GO:0023052,GO:0030097,GO:0030154,GO:0032501,GO:0032502,GO:0034754,GO:0035239,GO:0035295,GO:0035326,GO:0036211,GO:0042176,GO:0042221,GO:0042445,GO:0043170,GO:0043226,GO:0043227,GO:0043229,GO:0043231,GO:0043412,GO:0044212,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044424,GO:0044464,GO:0045137,GO:0045732,GO:0046545,GO:0046660,GO:0048008,GO:0048513,GO:0048514,GO:0048518,GO:0048519,GO:0048534,GO:0048608,GO:0048705,GO:0048731,GO:0048745,GO:0048856,GO:0048869,GO:0050789,GO:0050794,GO:0050896,GO:0051171,GO:0051173,GO:0051246,GO:0051247,GO:0051716,GO:0060021,GO:0060255,GO:0060322,GO:0060323,GO:0060324,GO:0060325,GO:0060537,GO:0061458,GO:0065007,GO:0065008,GO:0070887,GO:0071310,GO:0071407,GO:0071704,GO:0072001,GO:0072358,GO:0072359,GO:0080090,GO:0097159,GO:1901360,GO:1901363,GO:1901564
## 6 GO:0000003,GO:0003674,GO:0005215,GO:0005452,GO:0005488,GO:0005515,GO:0005575,GO:0005623,GO:0005886,GO:0005887,GO:0006810,GO:0006811,GO:0006812,GO:0006814,GO:0006820,GO:0006821,GO:0006873,GO:0006885,GO:0007275,GO:0007276,GO:0007283,GO:0008150,GO:0008324,GO:0008509,GO:0008510,GO:0008514,GO:0009987,GO:0015075,GO:0015077,GO:0015081,GO:0015103,GO:0015106,GO:0015108,GO:0015291,GO:0015293,GO:0015294,GO:0015297,GO:0015301,GO:0015318,GO:0015370,GO:0015672,GO:0015698,GO:0015701,GO:0015711,GO:0016020,GO:0016021,GO:0016323,GO:0016324,GO:0019725,GO:0019899,GO:0019953,GO:0022414,GO:0022804,GO:0022857,GO:0022890,GO:0030001,GO:0030003,GO:0030004,GO:0030641,GO:0031224,GO:0031226,GO:0032501,GO:0032502,GO:0032504,GO:0034220,GO:0035295,GO:0035725,GO:0042592,GO:0044425,GO:0044459,GO:0044464,GO:0044703,GO:0045177,GO:0046873,GO:0048232,GO:0048565,GO:0048609,GO:0048731,GO:0048856,GO:0048878,GO:0050801,GO:0051179,GO:0051234,GO:0051453,GO:0051704,GO:0055067,GO:0055080,GO:0055082,GO:0055085,GO:0055123,GO:0065007,GO:0065008,GO:0071702,GO:0071944,GO:0098590,GO:0098655,GO:0098656,GO:0098660,GO:0098661,GO:0098662,GO:0098771,GO:0099516,GO:1902476
## GO.description GS.Flat GS.Slope
## 1 thioredoxin-disulfide reductase activity 0.57190344 -0.57190344
## 2 - -0.29615673 0.29615673
## 3 carbonate dehydratase activity 0.82284313 -0.82284313
## 4 belongs to the actin family -0.44984287 0.44984287
## 5 TCDD-inducible poly ADP-ribose polymerase 0.04277864 -0.04277864
## 6 sodium:bicarbonate symporter activity 0.42124084 -0.42124084
## p.GS.Flat p.GS.Slope A.blue p.A.blue A.black p.A.black
## 1 3.296265e-05 3.296265e-05 0.6816236 1.839243e-07 -0.5833475 2.093028e-05
## 2 4.566951e-02 4.566951e-02 -0.4736549 8.846259e-04 0.4149787 4.135812e-03
## 3 2.282405e-12 2.282405e-12 0.9065745 4.306626e-18 -0.8069304 1.271757e-11
## 4 1.709469e-03 1.709469e-03 -0.5188606 2.204480e-04 0.5872458 1.785967e-05
## 5 7.777263e-01 7.777263e-01 -0.1544002 3.055833e-01 0.2095712 1.621606e-01
## 6 3.552683e-03 3.552683e-03 0.5190538 2.190498e-04 -0.6416509 1.544375e-06
## A.pink p.A.pink A.red p.A.red A.salmon p.A.salmon
## 1 0.27842751 6.097779e-02 -0.4083356 4.844425e-03 -0.7472873 2.437181e-09
## 2 -0.34099292 2.039159e-02 0.2228543 1.365748e-01 0.1622026 2.814859e-01
## 3 0.63501061 2.135745e-06 -0.6408045 1.610223e-06 -0.7482656 2.262822e-09
## 4 -0.39257347 6.963855e-03 0.1921334 2.008242e-01 0.3122323 3.464228e-02
## 5 -0.06841708 6.514225e-01 0.1073470 4.776535e-01 0.2320335 1.207387e-01
## 6 0.75371159 1.487555e-09 -0.6357069 2.065103e-06 -0.4624218 1.214176e-03
## A.greenyellow p.A.greenyellow A.turquoise p.A.turquoise A.cyan
## 1 0.6739431 2.838963e-07 -0.3330385 2.372177e-02 0.4267480
## 2 -0.1594412 2.898678e-01 0.5857920 1.895275e-05 -0.1271445
## 3 0.6280105 2.981403e-06 -0.4792416 7.526735e-04 0.5051039
## 4 -0.2841205 5.566998e-02 0.7612393 8.180878e-10 -0.4620600
## 5 -0.1447235 3.372496e-01 0.2102914 1.606897e-01 -0.1723410
## 6 0.0948033 5.308520e-01 -0.1847180 2.190911e-01 0.2164040
## p.A.cyan A.brown p.A.brown A.lightcyan p.A.lightcyan A.tan
## 1 0.0031008521 0.1271912 3.995972e-01 0.09110039 0.54709232 0.2176203
## 2 0.3997716440 -0.5933724 1.386015e-05 0.18659720 0.21435658 0.1819197
## 3 0.0003434504 0.3135143 3.386687e-02 0.22821859 0.12714273 0.1574452
## 4 0.0012264089 -0.7407478 3.969181e-09 0.33357821 0.02348224 0.0636769
## 5 0.2520848072 -0.3626069 1.326528e-02 0.18378842 0.22145957 0.4510512
## 6 0.1485953064 0.2471506 9.773671e-02 0.32845822 0.02583794 -0.2762368
## p.A.tan A.green p.A.green A.purple p.A.purple A.magenta
## 1 0.146270907 0.09661801 0.5229802445 0.08334567 5.818526e-01 -0.33576335
## 2 0.226274419 -0.21911333 0.1434547172 -0.15787549 2.946916e-01 -0.19726405
## 3 0.296026422 0.10379593 0.4924217545 0.29166711 4.921348e-02 -0.08868698
## 4 0.674179118 -0.50372335 0.0003587078 -0.06971277 6.452552e-01 -0.13140604
## 5 0.001655126 0.09486416 0.5305870974 -0.22683689 1.295240e-01 0.16146684
## 6 0.063125044 -0.23054378 0.1232098351 0.80793923 1.145965e-11 0.12353209
## p.A.magenta A.midnightblue p.A.midnightblue A.grey60 p.A.grey60
## 1 0.02253312 0.009700231 0.94898532 -0.0273657 0.856737263
## 2 0.18883104 0.100592699 0.50594353 -0.1887542 0.209010743
## 3 0.55780330 0.045801753 0.76245981 0.1449134 0.336609140
## 4 0.38402987 -0.094778219 0.53096122 -0.3901718 0.007348835
## 5 0.28370346 0.338323606 0.02146223 -0.2737761 0.065608424
## 6 0.41340438 -0.124918253 0.40814208 0.2161768 0.149032326glmmseq_exp <-read.csv("../../../output/Slope_Base/model_expression_prediction_allgenes.csv")
glmmseq_exp <- plyr::rename(glmmseq_exp, c("X"="Pocillopora_acuta_best_hit"))biomin_exp <- merge(biomin_mod, glmmseq_exp, by=c("Pocillopora_acuta_best_hit"), all.x=T)write.csv(biomin_exp, "../../../output/WGCNA/Biomineralization_toolkit_glmmseq_expression.csv")biomin_all <-read.csv("../../../output/Biomin_blast_Pocillopora_acuta_best_hit_glmmSeq.csv")
head(biomin_all)## Gene Dispersion AIC logLik
## 1 Pocillopora_acuta_HIv2___RNAseq.g13823.t1 0.29270610 366.4332 -177.2166
## 2 Pocillopora_acuta_HIv2___RNAseq.g13823.t1 0.29270610 366.4332 -177.2166
## 3 Pocillopora_acuta_HIv2___RNAseq.g13823.t1 0.29270610 366.4332 -177.2166
## 4 Pocillopora_acuta_HIv2___RNAseq.g25351.t1 0.13435721 815.9720 -401.9860
## 5 Pocillopora_acuta_HIv2___RNAseq.g7085.t1 0.25162583 491.9378 -239.9689
## 6 Pocillopora_acuta_HIv2___RNAseq.g22851.t1 0.02235176 759.7088 -373.8544
## meanExp X.Intercept. TreatmentVariable OriginFlat
## 1 4.316132 3.079614 0.16922532 0.07545055
## 2 4.316132 3.079614 0.16922532 0.07545055
## 3 4.316132 3.079614 0.16922532 0.07545055
## 4 12.237463 8.288576 -0.01709198 0.40440178
## 5 6.477549 4.543751 -0.19119978 0.09230384
## 6 11.722317 8.185091 -0.11783143 -0.07240262
## TreatmentVariable.OriginFlat se_.Intercept. se_TreatmentVariable
## 1 -0.3792527 0.1679989 0.2363277
## 2 -0.3792527 0.1679989 0.2363277
## 3 -0.3792527 0.1679989 0.2363277
## 4 0.2217089 0.1059042 0.1497752
## 5 0.3564220 0.1641296 0.2123867
## 6 0.1972383 0.0571238 0.0618745
## se_OriginFlat se_TreatmentVariable.OriginFlat Chisq_Treatment Chisq_Origin
## 1 0.24229985 0.34311487 0.003896091 0.4390691
## 2 0.24229985 0.34311487 0.003896091 0.4390691
## 3 0.24229985 0.34311487 0.003896091 0.4390691
## 4 0.15311285 0.21653465 0.676680451 22.6483669
## 5 0.22652466 0.30543197 0.013908932 2.6586173
## 6 0.08256803 0.08945199 0.275523829 0.1398413
## Chisq_Treatment.Origin Df_Treatment Df_Origin Df_Treatment.Origin P_Treatment
## 1 1.221739 1 1 1 0.9502294
## 2 1.221739 1 1 1 0.9502294
## 3 1.221739 1 1 1 0.9502294
## 4 1.048362 1 1 1 0.4107321
## 5 1.361758 1 1 1 0.9061183
## 6 4.861862 1 1 1 0.5996502
## P_Origin P_Treatment.Origin Treatment Origin Treatment.Origin
## 1 5.075721e-01 0.26901967 1 0.5306740470 0.9994279
## 2 5.075721e-01 0.26901967 1 0.5306740470 0.9994279
## 3 5.075721e-01 0.26901967 1 0.5306740470 0.9994279
## 4 1.945255e-06 0.30588455 1 0.0001683701 0.9994279
## 5 1.029902e-01 0.24323297 1 0.2466098058 0.9994279
## 6 7.084388e-01 0.02745669 1 0.6038776397 0.9994279
## Stable_OriginFC Variable_OriginFC maxGroupFC col RF13B
## 1 -0.1042361 0.4192811 Variable Not Significant 21.57376
## 2 -0.1042361 0.4192811 Variable Not Significant 21.57376
## 3 -0.1042361 0.4192811 Variable Not Significant 21.57376
## 4 -0.5833078 -0.9031151 Variable q_Origin < 0.05 5832.10758
## 5 -0.1318273 -0.6407294 Variable Not Significant 67.80326
## 6 0.1044247 -0.1800468 Variable Not Significant 4739.03686
## RF13D RF14B RF14C RF15B RF15D RF17B RF17D
## 1 25.50848 16.94826 19.20733 28.23482 21.29134 25.87585 14.01704
## 2 25.50848 16.94826 19.20733 28.23482 21.29134 25.87585 14.01704
## 3 25.50848 16.94826 19.20733 28.23482 21.29134 25.87585 14.01704
## 4 4927.06132 13021.91653 6113.26692 5054.03248 7484.79414 8321.89749 7240.96832
## 5 70.63887 54.61107 99.23788 164.36698 112.66669 121.50398 91.11075
## 6 4597.41325 2757.85926 3654.72843 3468.84911 3098.77752 2764.21551 2954.09080
## RF18B RF18D RF19B RF19C RF20B RF20C RF22B
## 1 22.30905 43.47845 0.0000 24.20314 34.69447 16.42775 13.92415
## 2 22.30905 43.47845 0.0000 24.20314 34.69447 16.42775 13.92415
## 3 22.30905 43.47845 0.0000 24.20314 34.69447 16.42775 13.92415
## 4 5891.93841 5089.20850 8016.9416 8866.41520 4982.12595 8081.53904 7060.31604
## 5 123.28687 101.44972 118.3772 140.60869 115.64823 135.98524 42.54600
## 6 2417.59689 3150.51549 3524.0588 2869.80032 3728.49908 2716.05423 4221.33720
## RF22C RF23A RF23C RF24B RF24D RF25A RF25C
## 1 18.19144 18.60133 36.2419 18.57332 23.35614 10.69681 16.36654
## 2 18.19144 18.60133 36.2419 18.57332 23.35614 10.69681 16.36654
## 3 18.19144 18.60133 36.2419 18.57332 23.35614 10.69681 16.36654
## 4 3512.96929 6530.04515 4346.0891 4820.26566 5124.51037 6260.69212 4758.11817
## 5 46.48924 249.64940 184.1480 184.75567 70.06842 92.45103 98.19927
## 6 3153.18302 3553.83267 3132.4753 3154.53200 3760.33872 3892.11199 3634.28212
## RS11B RS11D RS12A RS12C RS13A RS13C RS14B
## 1 37.03275 14.09808 0.00000 8.755265 20.55198 15.25056 29.16535
## 2 37.03275 14.09808 0.00000 8.755265 20.55198 15.25056 29.16535
## 3 37.03275 14.09808 0.00000 8.755265 20.55198 15.25056 29.16535
## 4 4408.84668 3804.31239 2168.53312 3681.588844 2701.11675 2417.21364 3177.56482
## 5 75.04005 300.39753 20.40972 12.257371 104.71721 80.06544 85.30865
## 6 3502.51878 3288.10581 3286.98596 2969.785817 3332.35599 3746.87177 3115.58845
## RS14C RS15B RS15D RS1B RS1C RS2B RS2C
## 1 14.95038 17.20235 15.77586 6.823476 33.59846 34.30448 35.76535
## 2 14.95038 17.20235 15.77586 6.823476 33.59846 34.30448 35.76535
## 3 14.95038 17.20235 15.77586 6.823476 33.59846 34.30448 35.76535
## 4 3123.75088 2933.85998 3892.69273 3421.973109 3717.90918 6782.61209 5170.97795
## 5 73.87249 85.15161 46.01292 55.725052 95.99559 59.81295 85.37536
## 6 3313.70871 3593.56992 3662.62815 3472.011931 3480.80008 2515.66212 3603.07098
## RS3B RS3D RS6A RS6D RS7B RS7C RS8B
## 1 44.72649 29.74233 29.32862 38.32429 30.78679 17.29086 9.894721
## 2 44.72649 29.74233 29.32862 38.32429 30.78679 17.29086 9.894721
## 3 44.72649 29.74233 29.32862 38.32429 30.78679 17.29086 9.894721
## 4 5164.93773 4277.13888 5058.41471 4024.05055 2932.44205 5029.82092 3791.477100
## 5 118.62244 77.71383 61.74446 120.29569 71.46934 82.81413 52.172164
## 6 2804.15671 2625.00044 3068.69956 3619.51637 3201.82649 3680.22360 3224.779455
## RS8C RS9A RS9C accessionnumber.geneID
## 1 14.95208 47.19477 19.44197 aug_v2a.09809.t1
## 2 14.95208 47.19477 19.44197 P13_g6918
## 3 14.95208 47.19477 19.44197 PFX18785.1
## 4 4200.60065 4097.88734 3391.45658 XP_022794351.1
## 5 75.69492 145.03759 102.65358 XP_022799541.1
## 6 3453.93103 2906.50718 4526.08973 P4_g9861
## definition
## 1 Mucin4-like protein
## 2 Sushi domain-containing
## 3 Mucin-4 [Stylophora pistillata]
## 4 mammalian ependymin-related protein 1-like [Stylophora pistillata]
## 5 uncharacterized protein LOC111337489 [Stylophora pistillata]
## 6 Viral inclusion protein
## Ref
## 1 Takeuchi et al., 2016
## 2 Drake et al., 2013
## 3 Peled et al., 2020
## 4 Peled et al., 2020
## 5 Peled et al., 2020
## 6 Drake et al., 2013library(tidyr)
biomin_all_filtered_long <- pivot_longer(biomin_all, cols=30:75, names_to = "Colony", values_to = "Expression")
biomin_all_filtered_long$Colony <- as.factor(biomin_all_filtered_long$Colony)
head(biomin_all_filtered_long)## # A tibble: 6 × 34
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 2 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 3 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 4 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 5 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 6 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## # ℹ 27 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <dbl>, Origin <dbl>,
## # Treatment.Origin <dbl>, Stable_OriginFC <dbl>, Variable_OriginFC <dbl>, …biomin_all_filtered_long <- biomin_all_filtered_long %>%
separate(Colony, into = c('Origin', 'Colony.number'), sep = 2)
library(stringr)
biomin_all_filtered_long$Colony <- as.numeric(str_extract(biomin_all_filtered_long$Colony.number, "[0-9]+"))
biomin_all_filtered_long <- biomin_all_filtered_long %>%
mutate(Treatment = trimws(str_remove(biomin_all_filtered_long$Colony.number, "(\\s+[A-Za-z]+)?[0-9-]+")))
biomin_all_filtered_long$Origin <- as.factor(biomin_all_filtered_long$Origin)
biomin_all_filtered_long$Treatment <- as.factor(biomin_all_filtered_long$Treatment)
biomin_all_filtered_long <- biomin_all_filtered_long %>%
mutate(Treatment2 = ifelse(Treatment == "A" | Treatment == "B", "Variable",
ifelse(Treatment == "C" | Treatment == "D", "Stable", NA)))
biomin_all_filtered_long$Treatment2 <- as.factor(biomin_all_filtered_long$Treatment2)
head(biomin_all_filtered_long)## # A tibble: 6 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 2 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 3 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 4 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 5 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## 6 Pocillopora_ac… 0.293 366. -177. 4.32 3.08 0.169
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>,
## # Treatment.Origin <dbl>, Stable_OriginFC <dbl>, Variable_OriginFC <dbl>, …#MEBlue signficant genes ##Thioredoxin reductase 1
biomin_all_filtered_long_g10093 <- biomin_all_filtered_long %>% filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g10093.t2")
biomin_all_filtered_long_g10093 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 2 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 3 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 4 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 5 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 6 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 7 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 8 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 9 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## 10 Pocillopora_a… 0.0876 385. -187. 5.74 3.88 -0.142
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:IRanges':
##
## collapse## The following object is masked from 'package:dplyr':
##
## collapselibrary(emmeans)## Warning: package 'emmeans' was built under R version 4.3.1g10093.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g10093 , na.action=na.exclude)
car::Anova(g10093.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 205.6032 1 < 2.2e-16 ***
## Origin 10.8768 1 0.0009738 ***
## Treatment2 0.2141 1 0.6435866
## Origin:Treatment2 1.3642 1 0.2428065
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3 <- emmeans(g10093.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 66.4 3.89 19 58.2 74.5
## RS 43.9 3.76 19 36.0 51.7
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 22.5 4.63 23 4.862 0.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g10093_sum<-summarySE(biomin_all_filtered_long_g10093 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g10093_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 63.79961 16.37758 4.938025 11.002606
## 2 RF Variable 11 65.98269 22.07239 6.655076 14.828433
## 3 RS Stable 12 47.40346 10.98308 3.170542 6.978315
## 4 RS Variable 12 41.95704 10.53729 3.041854 6.695076###Figure
pd<- position_dodge(0.2)
g10093_fig<-ggplot(data=g10093_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g10093,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Thioredoxin~reductase~1))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg10093_fig
##Carbonic anhydrase - STPCA2-2
biomin_all_filtered_long_g13824 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g13824.t1")
biomin_all_filtered_long_g13824 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 2 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 3 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 4 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 5 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 6 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 7 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 8 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 9 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## 10 Pocillopora_a… 2.35 453. -220. 4.66 2.36 -0.0648
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g13824.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g13824 , na.action=na.exclude)
car::Anova(g13824.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 77.9228 1 <2e-16 ***
## Origin 78.0753 1 <2e-16 ***
## Treatment2 2.0231 1 0.1549
## Origin:Treatment2 1.0392 1 0.3080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g13824.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 158.76 17.7 19 121.7 195.8
## RS -5.68 17.3 19 -41.9 30.5
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 164 17 23 9.671 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g13824_sum<-summarySE(biomin_all_filtered_long_g13824 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g13824_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 156.36191 109.769844 33.096853 73.744384
## 2 RF Variable 11 135.61447 87.446714 26.366176 58.747502
## 3 RS Stable 12 10.39836 9.661042 2.788902 6.138333
## 4 RS Variable 12 10.23764 8.081504 2.332929 5.134743###Figure
pd<- position_dodge(0.2)
g13824_fig<-ggplot(data=g13824_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g13824,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Carbonic~anhydrase~("STPCA2-2")), limits=c(0,300))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg13824_fig## Warning: Removed 1 rows containing missing values (`geom_point()`).
##Mammalian ependymin-related protein 1-like
biomin_all_filtered_long_g25351 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g25351.t1")
biomin_all_filtered_long_g25351 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 2 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 3 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 4 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 5 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 6 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 7 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 8 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 9 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## 10 Pocillopora_a… 0.134 816. -402. 12.2 8.29 -0.0171
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g25351.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g25351 , na.action=na.exclude)
car::Anova(g25351.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 151.907 1 < 2.2e-16 ***
## Origin 10.058 1 0.001517 **
## Treatment2 1.959 1 0.161619
## Origin:Treatment2 1.039 1 0.308064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g25351.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 6454 354 19 5713 7195
## RS 3869 339 19 3159 4578
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2585 482 23 5.359 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g25351_sum<-summarySE(biomin_all_filtered_long_g25351 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g25351_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 5958.631 1717.4980 517.8451 1153.8309
## 2 RF Variable 11 6890.207 2342.7581 706.3681 1573.8863
## 3 RS Stable 12 3894.293 756.1492 218.2815 480.4342
## 4 RS Variable 12 3886.639 1300.8557 375.5247 826.5242###Figure
pd<- position_dodge(0.2)
g25351_fig<-ggplot(data=g25351_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g25351,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Mammalian~ependymin-related~protein~1-like),limits=c(2000,10000))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg25351_fig## Warning: Removed 1 rows containing missing values (`geom_point()`).
##Cephalotoxin-like protein
biomin_all_filtered_long_g5013 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g5013.t1")
biomin_all_filtered_long_g5013 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 2 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 3 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 4 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 5 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 6 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 7 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 8 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 9 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## 10 Pocillopora_a… 0.471 304. -146. 3.16 1.83 0.206
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g5013.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g5013 , na.action=na.exclude)
car::Anova(g5013.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 45.2530 1 1.732e-11 ***
## Origin 6.6236 1 0.01006 *
## Treatment2 0.0840 1 0.77190
## Origin:Treatment2 0.0561 1 0.81279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g5013.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 13.55 1.49 19 10.43 16.7
## RS 7.14 1.43 19 4.16 10.1
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 6.41 1.98 23 3.237 0.0036
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g5013_sum<-summarySE(biomin_all_filtered_long_g5013 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g5013_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 13.339697 7.878658 2.375505 5.292955
## 2 RF Variable 11 14.078906 6.431123 1.939056 4.320487
## 3 RS Stable 12 6.188827 5.219635 1.506779 3.316398
## 4 RS Variable 12 7.764125 6.413793 1.851503 4.075130###Figure
pd<- position_dodge(0.2)
g5013_fig<-ggplot(data=g5013_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g5013,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Cephalotoxin-like~protein), breaks=c(0,5,10,15,20,25))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg5013_fig
#MEPink signficant genes ##Carbonic anhydrase - STPCA2-1
biomin_all_filtered_long_g12304 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g12304.t1")
biomin_all_filtered_long_g12304 ## # A tibble: 184 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 2 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 3 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 4 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 5 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 6 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 7 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 8 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 9 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## 10 Pocillopora_a… 0.457 841. -415. 11.6 7.88 0.0448
## # ℹ 174 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g12304.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g12304 , na.action=na.exclude)
car::Anova(g12304.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 196.964 1 < 2.2e-16 ***
## Origin 55.461 1 9.531e-14 ***
## Treatment2 14.008 1 0.0001820 ***
## Origin:Treatment2 11.399 1 0.0007348 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g12304.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 5238 384 19 4434 6042
## RS 2770 375 19 1985 3555
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2468 371 161 6.657 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g12304_sum<-summarySE(biomin_all_filtered_long_g12304 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g12304_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 44 5960.191 2423.309 365.3276 736.7533
## 2 RF Variable 44 4766.484 2433.696 366.8934 739.9111
## 3 RS Stable 48 2494.782 1512.223 218.2706 439.1038
## 4 RS Variable 48 2791.885 1392.250 200.9540 404.2674###Figure
pd<- position_dodge(0.2)
g12304_fig<-ggplot(data=g12304_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g12304,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Carbonic~anhydrase~("STPCA2-1")), limits=c(0,10000))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg12304_fig## Warning: Removed 4 rows containing missing values (`geom_point()`).
##SLC4 gamma- solute carrier family 4 member gamma
biomin_all_filtered_long_g15280 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g15280.t1")
biomin_all_filtered_long_g15280 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 2 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 3 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 4 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 5 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 6 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 7 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 8 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 9 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## 10 Pocillopora_a… 0.376 425. -206. 5.12 3.45 0.0890
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g15280.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g15280 , na.action=na.exclude)
car::Anova(g15280.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 82.1177 1 < 2.2e-16 ***
## Origin 9.0855 1 0.002576 **
## Treatment2 0.7652 1 0.381698
## Origin:Treatment2 0.9266 1 0.335743
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g15280.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 52.0 4.36 19 42.9 61.2
## RS 32.1 4.17 19 23.4 40.9
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 19.9 6.03 23 3.300 0.0031
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g15280_sum<-summarySE(biomin_all_filtered_long_g15280 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g15280_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 55.84385 24.85863 7.495159 16.700254
## 2 RF Variable 11 48.22013 24.68527 7.442889 16.583790
## 3 RS Stable 12 30.12777 14.24152 4.111172 9.048629
## 4 RS Variable 12 34.11841 16.62673 4.799723 10.564118###Figure
pd<- position_dodge(0.2)
g15280_fig<-ggplot(data=g15280_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g15280,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Solute~carrier~family~4~member~gamma))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg15280_fig
##SLC4A7 - sodium bicarbonate cotransporter 3-like isoform X2
biomin_all_filtered_long_g7402 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g7402.t1")
biomin_all_filtered_long_g7402 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 2 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 3 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 4 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 5 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 6 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 7 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 8 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 9 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## 10 Pocillopora_a… 0.229 640. -314. 8.97 6.14 0.0129
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g7402.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g7402 , na.action=na.exclude)
car::Anova(g7402.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 87.9465 1 <2e-16 ***
## Origin 6.3555 1 0.0117 *
## Treatment2 0.8936 1 0.3445
## Origin:Treatment2 0.8273 1 0.3630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g7402.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 660 54.9 19 545 775
## RS 464 52.6 19 353 574
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 197 75 23 2.621 0.0153
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g7402_sum<-summarySE(biomin_all_filtered_long_g7402 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g7402_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 710.2028 264.1958 79.65802 177.4891
## 2 RF Variable 11 611.9007 312.5396 94.23425 209.9670
## 3 RS Stable 12 445.6090 228.0413 65.82984 144.8905
## 4 RS Variable 12 478.2523 190.4959 54.99144 121.0353###Figure
pd<- position_dodge(0.2)
g7402_fig<-ggplot(data=g7402_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g7402,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(sodium~bicarbonate~cotransporter~3-like~isoform~X2))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg7402_fig
##Protein lingerer-like
biomin_all_filtered_long_g7908 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g7908.t1")
biomin_all_filtered_long_g7908 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 2 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 3 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 4 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 5 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 6 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 7 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 8 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 9 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## 10 Pocillopora_a… 0.138 426. -207. 5.95 3.98 -0.0743
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g7908.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g7908 , na.action=na.exclude)
car::Anova(g7908.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 207.4212 1 < 2.2e-16 ***
## Origin 15.8613 1 6.816e-05 ***
## Treatment2 0.1605 1 0.6887
## Origin:Treatment2 0.0058 1 0.9392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g7908.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 81.2 4.07 19 72.7 89.7
## RS 49.9 3.89 19 41.8 58.0
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 31.3 5.63 23 5.556 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g7908_sum<-summarySE(biomin_all_filtered_long_g7908 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g7908_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 82.80022 22.33628 6.734642 15.005717
## 2 RF Variable 11 79.54283 21.44921 6.467180 14.409774
## 3 RS Stable 12 51.10111 16.98594 4.903419 10.792352
## 4 RS Variable 12 48.70220 15.09645 4.357969 9.591824###Figure
pd<- position_dodge(0.2)
g7908_fig<-ggplot(data=g7908_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g7908,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(protein~lingerer-like))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg7908_fig
#Signficant genes not in WGCNA modules ##Vitellogenin
biomin_all_filtered_long_g13222 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g13222.t1b")
biomin_all_filtered_long_g13222 ## # A tibble: 184 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 2 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 3 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 4 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 5 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 6 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 7 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 8 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 9 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## 10 Pocillopora_a… 0.0506 621. -305. 9.50 6.53 -0.123
## # ℹ 174 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g13222.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g13222 , na.action=na.exclude)
car::Anova(g13222.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 1126.7286 1 < 2.2e-16 ***
## Origin 41.2035 1 1.372e-10 ***
## Treatment2 2.9349 1 0.08669 .
## Origin:Treatment2 0.5477 1 0.45926
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g13222.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 832 22.8 19 784 879
## RS 636 22.2 19 590 683
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 195 24.2 161 8.053 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g13222_sum<-summarySE(biomin_all_filtered_long_g13222 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g13222_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 44 859.5886 158.58611 23.90776 48.21459
## 2 RF Variable 44 821.8027 114.92648 17.32582 34.94084
## 3 RS Stable 48 660.1489 146.62900 21.16407 42.57662
## 4 RS Variable 48 599.7645 70.60441 10.19087 20.50138###Figure
pd<- position_dodge(0.2)
g13222_fig<-ggplot(data=g13222_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g13222,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Vitellogenin))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg13222_fig
##uncharacterized protein LOC111323869
biomin_all_filtered_long_g21232 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g21232.t1")
biomin_all_filtered_long_g21232 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 2 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 3 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 4 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 5 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 6 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 7 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 8 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 9 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## 10 Pocillopora_a… 0.0605 353. -171. 5.22 3.54 -0.0600
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g21232.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g21232 , na.action=na.exclude)
car::Anova(g21232.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 209.9697 1 < 2e-16 ***
## Origin 4.7421 1 0.02943 *
## Treatment2 0.1307 1 0.71769
## Origin:Treatment2 0.0003 1 0.98700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g21232.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 42.0 2.4 19 36.9 47.0
## RS 33.7 2.3 19 28.8 38.5
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 8.29 3.02 23 2.747 0.0115
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g21232_sum<-summarySE(biomin_all_filtered_long_g21232 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g21232_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 42.29276 10.185331 3.070993 6.842599
## 2 RF Variable 11 41.06536 11.600467 3.497672 7.793300
## 3 RS Stable 12 33.84473 7.075632 2.042559 4.495642
## 4 RS Variable 12 32.69394 10.921097 3.152649 6.938934###Figure
pd<- position_dodge(0.2)
g21232_fig<-ggplot(data=g21232_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g21232,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(uncharacterized~protein~LOC111323869))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg21232_fig
##uncharacterized protein LOC111345150
biomin_all_filtered_long_g20587 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g20587.t2")
biomin_all_filtered_long_g20587 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 2 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 3 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 4 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 5 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 6 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 7 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 8 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 9 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## 10 Pocillopora_a… 4.85 271. -129. 1.93 2.59 0.106
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g20587.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g20587 , na.action=na.exclude)
car::Anova(g20587.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 0.1659 1 0.683767
## Origin 7.2048 1 0.007271 **
## Treatment2 0.1684 1 0.681502
## Origin:Treatment2 0.0040 1 0.949524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g20587.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 2.06 2.78 19 -3.77 7.88
## RS 13.31 2.68 19 7.70 18.93
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -11.3 3.35 23 -3.355 0.0027
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g20587_sum<-summarySE(biomin_all_filtered_long_g20587 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g20587_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 1.069422 2.379314 0.7173901 1.598445
## 2 RF Variable 11 2.503911 3.677310 1.1087507 2.470451
## 3 RS Stable 12 12.756815 15.369960 4.4369253 9.765607
## 4 RS Variable 12 14.497631 15.007069 4.3321675 9.535036###Figure
pd<- position_dodge(0.2)
g20587_fig<-ggplot(data=g20587_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g20587,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(uncharacterized~protein~LOC111345150), limits=c(0,30))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg20587_fig## Warning: Removed 4 rows containing missing values (`geom_point()`).
##Late embryogenesis protein
biomin_all_filtered_long_g16715 <- biomin_all_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g16715.t1")
biomin_all_filtered_long_g16715 ## # A tibble: 46 × 36
## Gene Dispersion AIC logLik meanExp X.Intercept. TreatmentVariable
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 2 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 3 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 4 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 5 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 6 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 7 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 8 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 9 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## 10 Pocillopora_a… 0.154 847. -417. 12.6 8.39 0.255
## # ℹ 36 more rows
## # ℹ 29 more variables: OriginFlat <dbl>, TreatmentVariable.OriginFlat <dbl>,
## # se_.Intercept. <dbl>, se_TreatmentVariable <dbl>, se_OriginFlat <dbl>,
## # se_TreatmentVariable.OriginFlat <dbl>, Chisq_Treatment <dbl>,
## # Chisq_Origin <dbl>, Chisq_Treatment.Origin <dbl>, Df_Treatment <int>,
## # Df_Origin <int>, Df_Treatment.Origin <int>, P_Treatment <dbl>,
## # P_Origin <dbl>, P_Treatment.Origin <dbl>, Treatment <fct>, …library(nlme)
g16715.lme <- lme(Expression~Origin*Treatment2, random = ~1|Colony, data=biomin_all_filtered_long_g16715 , na.action=na.exclude)
car::Anova(g16715.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Expression
## Chisq Df Pr(>Chisq)
## (Intercept) 265.1538 1 < 2.2e-16 ***
## Origin 32.7124 1 1.069e-08 ***
## Treatment2 0.3604 1 0.5483
## Origin:Treatment2 1.3445 1 0.2462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g16715.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 8304 381 19 7507 9101
## RS 4973 365 19 4209 5737
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 3331 505 23 6.602 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g16715_sum<-summarySE(biomin_all_filtered_long_g16715 , measurevar='Expression', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g16715_sum## Origin Treatment2 N Expression sd se ci
## 1 RF Stable 11 8076.757 1731.637 522.1082 1163.3295
## 2 RF Variable 11 8462.423 1729.677 521.5172 1162.0128
## 3 RS Stable 12 4230.613 1879.259 542.4954 1194.0243
## 4 RS Variable 12 5647.539 1237.158 357.1369 786.0529###Figure
pd<- position_dodge(0.2)
g16715_fig<-ggplot(data=g16715_sum, aes(y=Expression, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=biomin_all_filtered_long_g16715,aes(y=Expression, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Expression-se, ymax=Expression+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_color_manual("Origin", values=c("RS"='#26519B', "RF"= "#FE180C"))+
scale_y_continuous(expression(Late~embryogenesis~protein))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg16715_fig
#Comparison of all sig. genes
biomin_compare_figs<-cowplot::plot_grid(g13824_fig,g25351_fig,g10093_fig,g5013_fig, g12304_fig,g7908_fig,g12304_fig,g13222_fig,g16715_fig,g20587_fig,g21232_fig,g15280_fig,g7402_fig, nrow=4)## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).
## Removed 4 rows containing missing values (`geom_point()`).
## Removed 4 rows containing missing values (`geom_point()`).biomin_compare_figs